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Abstract 

A two-fold is a singular point on the discontinuity surface of a piecewise-smooth vector 
field, at which the vector field is tangent to the discontinuity surface on both sides. If an orbit 
passes through an invisible two-fold (also known as a Teixeira singularity) before settling to 
regular periodic motion, then the phase of that motion cannot be determined from initial 
conditions, and in the presence of small noise the asymptotic phase of a large number of 
sample solutions is highly random. In this paper we show how the probability distribution 
of the asymptotic phase depends on the global nonlinear dynamics. We also show how the 
phase of a smooth oscillator can be randomised by applying a simple discontinuous control 
law that generates an invisible two-fold. We propose that such a control law can be used 
to desynchronise a collection of oscillators, and that this manner of phase randomisation is 
fast compared to existing methods (which use fixed points as phase singularities) because 
there is no slowing of the dynamics near a two-fold. 


1 Introduction 

When a solution of a system of ordinary differential equations approaches a stable periodic or¬ 
bit, rather than its location in phase space at a given time, a more useful quantity is often 
the asymptotic (or eventual) phase of the orbit. The asymptotic phase distinguishes the long¬ 
term behaviour of different orbits, and sets of points with the same asymptotic phase are called 
isochrons. If different isochrons intersect they create a phaseless point, where the asymptotic 
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phase is undefined. Near a phaseless point the asymptotic phase is highly sensitive to pertur¬ 
bations. In smooth systems, phaseless points are often unstable equilibria which can only be 
approached asymptotically (in backward time). In piecewise-smooth systems there exist phase¬ 
less points that can be intersected in finite time, without any slowing of the dynamics suffered 
near an equilibrium. The two-fold singularity, an enigma of piecewise-smooth dynamical systems 
theory, is one such phaseless point. 

There are many reasons why it might be desirable to alter the phase of an existing oscillatory 
motion. One approach to tackling Parkinson’s disease, for instance, is to disturb the synchro¬ 
nised neural activity associated with physical tremors HIE]. Desynchronisation can be achieved 
with pulses [3], pulse trains [1], a control that utilises time-delay [5], or some other well-chosen 
feedback law [6]. Phase randomisation for prototypical neuron models is achieved in [7] by briefly 
applying a control that transports the orbit to the close proximity of a phaseless point. The orbit 
subsequently returns to periodic motion but now has a different asymptotic phase. This phase 
is highly sensitive to the precise point at which the orbit is located when the control is removed, 
so small stochasticity in the system causes the resulting asymptotic phases of different neurons 
to be highly randomised. 

In piecewise smooth dynamical systems, two-fold singularities were first described by Filippov 
[8], and have garnered interest as points where both the forward and backward time uniqueness 
of a flow can break down, in an otherwise deterministic flow [acniE]. For a vector field that 
is discontinuous on some hypersurface - the discontinuity surface - a two-fold is a point where 
the flow has quadratic tangencies (‘folds’) to both sides of the surface simultaneously. Two-folds 
occur generically at isolated points in three-dimensional piecewise-smooth vector fields, and on 
(n — 2)-dimensional manifolds in n-dimensional vector fields [11]. They have been identified 
in models of circuit systems [12] and mechanical systems [TB], and have deep connections to 
folded nodes of slow-fast systems [TT] . The dynamics local to a two-fold depends on whether the 
vector field on either side is curving toward or away from the discontinuity surface, and on the 
alignments of the two vector fields relative to each other m- 

Here we illustrate the practical implications of a two-fold for phase randomisation. As well 
as suggesting how two-folds might affect real systems, this suggests a use for them as control 
devices for the desynchronisation of oscillators. To this end we shall consider systems where 
a two-fold occurs either naturally in a piecewise-smooth system, or is introduced to a smooth 
system via a discontinuous control. We will simulate orbits in the presence of small noise as they 
pass close to a two-fold, and focus on how their phases are randomised by the nearby presence 
of the singularity. The two-fold is not a zero of the vector field, so there is no slowing of the flow 
during phase randomisation achieved in this way. 

The central result of the paper is summarised by Fig. [T] (to be produced in showing the 
phase distributions of 10^ sample solutions from a fixed initial point, simulated passing through 
a two-fold singularity. The central result of the paper is summarised by Fig. [1] (to be produced in 
Q, showing the phase distributions of 10"^ sample solutions from a fixed initial point, simulated 
passing through a two-fold singularity. To produce these figures we shall simulate the normal form 
of the two-fold, add higher order terms such that orbits emanating from the two-fold approach a 
stable periodic orbit, and add small noise. Fig. [T]-A shows a histogram of the phase (j)T-, computed 
relative to a reference time T and limiting to the asymptotic phase as T —)■ cxd. This shows an 
apparently uniform distribution, hence a fully randomised phase of solutions from the same initial 
point. Fig. [T]-B shows an analogous histogram using different higher order terms, showing how 
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Figure 1: Each panel shows a histogram of the phase of 10^ sample solutions from the same 
initial condition, of the general stochastic system (14. Ih . The two panels correspond to different 
choices for the higher order terms which affect the global dynamics of (14. ip . One sample solution 
corresponding to each panel is shown in Fig. [7] (see caption of this hgure for the parameter values 
used). The solid curves show a theoretical approximation for the probability density function of 
0^ as derived in 1 14.21 


nonlinearity of the flow can be used to influence the distribution. The black curves show the 
theoretical approximations we derive for the probability density functions of the phase. 

Building on this, we propose that a periodic orbit in a smooth system can be made to take an 
excursion through a two-fold, by applying a control action, before returning to smooth periodic 
motion with a randomised phase. We provide simulations of a simple model in which the effect is 
to desynchronise a collection synchronised oscillators. As discussed in ^ two-folds have potential 
advantages over equilibria as phaseless points in the context of phase randomisation. Although 
we provide a motivating toy example, calculating and optimising realisable control actions (as in 
e.g. 0) are beyond the scope of this paper. 

In ^ we briefly set out the normal form equations of the two-fold singularity, and review 
a few pertinent details of the local dynamics, including a novel dehnition of polar coordinates 
that naturally £t with the local dynamics (which should be of particular interest for dynamicists 
studying two-folds). We include higher order terms with the normal form in ^ and extend these 
definitions before simulating a flow through the two-fold subject to stochastic perturbations in 
^ In §3 we give an example of an application to desynchronising smooth oscillators using a 
discontinuous control, with some closing remarks in §21 


2 Deterministic dynamics of the normal form 

There exist three main kinds of two-fold singularity, depending on whether both folds are visible 
(curving away the discontinuity surface), or invisible (curving toward the discontinuity surface), 
or one is visible and one is invisible. In each case there exist local conditions under which the 
flow traverses the singularity in hnite time. In this paper we focus on the case where both folds 
are invisible, called an invisible two-fold or Teixeira singularity, in which the entire local flow 
traverses the singularity in hnite time. The added interest of the invisible two-fold is that the 
local how winds repeatedly around the singularity, giving oscillatory behaviour. 

In three dimensions where X = {x,y,z), the normal form of the invisible two-fold is the 
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piecewise-smooth system 

Fl{X) , X < 0 , 

Fr{X) , X > 0 , 

where 

Fl{X) = (z, y-, 1) , Fr{X) = (- 7 /, 1, , 

and V-,V+ e M are parameters, [I5l[l6]. The two-fold is located at the origin, X = (0,0,0). 
The discontinuity surface x = 0 consists of attracting and repelling sliding regions denoted A 
and R (where the flow is conhned to the surface and so slides along it), and two crossing regions 
(where the flow passes transversally though the surface), illustrated in Fig. |2l 
The system is solved across the discontinuity by forming the Filippov system for fl2.ip . Using 
Cj to denote the coordinate vector in a Filippov solution of fl2.ip is an absolutely continuous 
function ip(t) that satishes 

^ = Fiicfit)) , if el(p{t) < 0 , 

^ = FR{Lp{t)) , if > 0 , (2.3) 

^ e {(1 - q)FL{ip{t)) + qFR{ip{t)) | 0 < g < 1} , if ejipit) = 0 , 

for almost all values of t [8]. We let (p{t] Xo,to) denote a Filippov solution to fl2.ip with initial 

condition (p{to) = Xq. Each ip{t;Xo,tQ) is a concatenation of several smooth orbit segments, 
including segments of evolution outside the discontinuity surface x = 0, and segments of evolution 



( 2 . 1 ) 

( 2 . 2 ) 



Figure 2: A schematic of the two-fold of fl2.ip . The discontinuity surface x = 0 is made up 
of an attracting sliding region A {y^z > 0), a repelling sliding region R {y,z < 0), and two 
crossing sliding regions, C~^ {y < 0, z > 0) and C~ {y > 0, z < 0). Parts of two orbits and their 
intersections with the discontinuity surface are shown. One orbit starts in R and has subsequent 
intersections in C~^ and C~. A second orbit starts and returns to an invariant line ( G C~, see 
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on a; = 0 referred to as sliding motion. Solving the convex combination in fl2.3p for i; = 0 gives 
the system that sliding motion obeys on A and R as 
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Various details of the sliding flow can be derived from fl2.4p . see [151 US]- Throughout this paper 
we assume 

V < 0 , V+ < 0 , V-V+ > 1 , (2.5) 

in which case a typical orbit of fl2.3p and fl2.4p will intersect the two-fold in forward time or back¬ 
ward time or both (see e.g. [I5]), giving the geometry in Fig. |2l With (I2.5p . the matrix in (12.4p 
has distinct negative eigenvalues, the smaller of which is A = -|- V~ + [(V’*' — -|- 4]^/^). 

As orbits of fl2.ip slide into the singularity inside A, or slide out of the singularity in R, they do 
so tangent to the weak eigenvector of the matrix (“weak” meaning associated with the smallest 
eigenvalue) and so are tangent to the line z = {X — V~)y on x = 0. The time to reach or depart 
the singularity is finite. Specifically, the time is t = (V“ — A — l)y/X from an initial point on the 
line 2 ; = (A — V~)y with x = 0. 

In the remainder of this section we review crossing solutions to (12.ip in §2.H introduce polar 
coordinates for describing the flow as it spirals away from the two-fold in 1 12.21 and define the 
notion of “viable” Filippov solutions in the context of simulation in 1 12.31 In later sections we 
apply these concepts to more general piecewise-smooth systems. 

2.1 Crossing dynamics and an unstable cone 

The left and right half systems of fl2.ip . X = Fl{X) and X = Fr(X), have solutions 

-^0; to) = + Zo{t — to) + -{t — to)^5 2/0 + 1 ^ {t — to), Zq + t — , (2.6) 

and 

^PR{t;Xo,to) = (^xo-yo{t-to)yo + t-to, zo + V^{t-to)'^ , (2.7) 

respectively, where Xq = (xo,yo,Zo)- In what follows we use fl2.6|) and fl2.7p to study Filippov 
solutions (p(t;0,y,z,to) of fl2.ip . Orbits wind around the regions A and R, making repeating 
crossings of x = 0 via the regions C^, which can be understood with a return map that has 
been studied in considerable detail usmi]. Here we review a few details needed for the ensuing 
analysis. 

First consider the region z < 0 on x = 0, where the flow cpi is directed away from the 
discontinuity surface. If ?/ > 0, then (0,y,z) G C~, see Fig. [2l (pR is directed toward the surface 
so (p(t;0,y,z,to) immediately enters x < 0. If instead y < 0, then (0,y,z) G R where both ipi 
and ipR are directed away from the surface, so ip{t; 0, y, z, R) is not uniquely determined (the may 
slide along x = 0 following fl2.4p . or enter either x < 0 or x > 0 at any instant). For definiteness, 
when ?/ < 0 we ignore solutions that follow fl2.4|) through R, and consider the Filippov solution 
ip(t;0,y, z,to) that immediately enters x < 0. From fl2.6l) we find that (p(t;0,y, z,to) resides in 
X < 0 until the later time t = R — 2z when it is located at (0, y — 2V~z, —z). 
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Second, consider y <0 and assnme ip{t;0,y, z,to) immediately enters a; > 0. From fl2.7p and 
applying analogous arguments to the case 2 ; < 0 , we find that ip{t; 0 , y, z, to) resides in a; > 0 until 
t = to — 2y when it is located at (0, —y, —2V^y + z). 

From these observations we can construct a return map capturing crossing dynamics, formu¬ 
lated in the following proposition (see [T 6 j for a complete proof). 


Proposition 2.1. Suppose V < 0, V > 1. For any ?/ G M and z < 0, let {0,y',z') 
denote the next intersection of(p{t; 0, y, z, to) with the discontinuity surface at a point with y' > 0. 
If z < y, then z' <t) and 


y' 


1 

1 

to 

1 


y 

z' 


-2V+ 4V-V+-1 


z 


the smooth orbit from {y, z) to {y', z') taking a time 

t-to = 2 ((2V~ - l)z - y) , 


( 2 . 8 ) 


(2.9) 


whereas if z > ^y_y+_^ y, then z' > 0. 

The return map (12.81) is linear and has a saddle-type fixed point at {y,z) = (0,0). The 
unstable multiplier (i.e. the eigenvalue of the matrix in (12.81) with modulus greater than 1 ) is 


/i 


2V-V+ 



- 1 , 


and the corresponding eigenvector is ( 1 , 7 ) where 



( 2 . 10 ) 


( 2 . 11 ) 


In the full three dimensional space occupied by fl2.ip . this eigenvector corresponds to a ray 


C = {(o,y,7y) I 2 / > 0 } , 


( 2 . 12 ) 


that emanates from the two-fold. 

By evolving the flow forward from ( according to 
is given implicitly by 


X = 




fl2.6l) and fl2.7p we obtain a surface A that 


a; < 0 
a; > 0 ’ 


(2.13) 


where 

, V+y^ -2V+V-yz + V-z^ 

=(V. =- 2(V^V--1) - 


(2.14) 


Fig. |3] shows a plot of A. A is non-differentiable at x = 0 and consists of one conical portion on 
each side of x = 0 . 
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Figure 3: A plot of 02.131) using V = —0.5 and V~^ = —2.5. This is the unstable manifold A 
02.13p . One Filippov solution 'ipait) is shown in blue. 


For any initial point on the ray (, the next intersection of (p(f; 0, ?/, 'yy, to) with C~ is (0, yy, y'jy), 
shown in Fig. [2l By substituting z = 'yy into 02.91) . we find that this iteration of 02.81) corresponds 
to an evolution time of 

f-to = 2 ((21/--1)7-1) I/. (2.15) 


By performing this iteration backward through repeated crossings, we can construct a Filippov 
solution that approaches the two-fold in backward time. By taking the limit of infinitely many 
intersections, we obtain a Filippov solution that emanates from the two-fold and has infinitely 
many intersections with C, (in finite time). Given any a > 0, let (defined for f > 0) denote 

the unique Filippov solution that is located at the two-fold at f = 0, and located on ^ at f = a. 
That is, '0a(O) = (0, 0, 0) and '0a(a) £ C- Such an orbit is shown in Fig. [3l 

By using fl2.15p and the classical formula for the sum of a geometric series, we determine the 
?/-value of '0a(o) fo be 


elV’a(a) = 


1 + 


1 -^ 


1 - 


V-V+ 


( 2 . 16 ) 


After t = a, 'ipait) next intersects ( at t = ya. Consequently ipa G (, for all n G Z. It 

follows that each is the same orbit as 'ipait), and so we can characterise all the 'ipait) by 

restricting our attention to values of a in an interval a* < a < ya*, for some a* > 0. 

We can therefore write 

^ = {i’ait) \ t > 0,1 < a < y} . (2-17) 


The intersection of A with C is Since ( corresponds to an unstable eigenvector of the crossing 
map, we refer to A as an unstable manifold of the two-fold. 


2.2 Polar coordinates 

On A, orbits of fl2.1l) spiral out from the two-fold, as in Fig. [3l For this reason it is natural 
to introduce polar coordinates, but the obvious spherical or cylindrical coordinates bear no 
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relation to the system’s behaviour and yield complicated expression that give no insight into the 
dynamical system. Instead let us seek radial and angular coordinates r(a;, y) and 6{x, y) that may 
be interpreted as “cylindrical polar coordinates” for orbits winding around the sliding regions, 
derived from the geometry of A such that the dynamical system takes a particularly simple form. 

To obtain a suitable notion of the phase 6*, we begin by finding the times > 0 and < 0 
that an orbit takes to travel from to a given point in a: < 0 and in x > 0, respectively. 

First, choose any x < 0 and ?/ G M. We let ?/o > 0 be such that by evolving fl2.ip forward 
from (0, j/o, 7?/o) we arrive at (x, ?/, z) G A, for some 2 ; G M, without again intersecting x = 0, and 
let tl > 0 he the corresponding evolution time. Note that the value of 2 ; is given from (I2.13p by 

X = ^S(y,z). 

We have 93 ^(tl; 0, yo, 72 /o, 0) = {x,y,z), and so by ([2J]), 

X = -yyaTL + ^'^1 , y = yo + V~tl , ^ = 7l/o + • (2.18) 

By solving the first two of these equations simultaneously for y^ and we obtain 


tl = 



(2.19) 


and yo = y - V~tl. 

Second, choose any x > 0 and y G M. We let yo > 0 be such that by evolving fl2.ip backward 
from (0, yo, ■yyo) we arrive at (x, y, z) G A, for some 2 ; G M, without hrst reintersecting x = 0, and 
let Ti? < 0 be the corresponding evolution time. Here 2 ; is given by x = .^S(?/,^), and in the 
same manner as above we obtain 

TR = y- ^y^ + 2 x , ( 2 . 20 ) 

and yo = y- tr. 

Using these we can now introduce a set of polar coordinates. We let the positive y-axis 
correspond to 0 = 0 and r = y, i.e. 

r( 0 , y) = y, 0 ( 0 , |/) = 0 , for all |/ > 0 . ( 2 . 21 ) 

The positive j/-axis corresponds to points on (. Naturally we want to dehne 6 such that a change 
of 271 corresponds to one complete revolution from ( back to itself. By fl2.15p . forward evolution 
from (0, y, 7 ?/) returns to C in a time proportional to y. This suggests that we want 9 to be 
inversely proportional to y (or the distance from the two-fold, r). Moreover, the orbit returns to 
C at the point {0, fiy, y.jy). That is, the value of r increases from y to y.y and so changes by an 
amount proportional to y, suggesting r should be constant. In summary, we would like to define 
r and 6 so that they satisfy fl 2 . 2 ip and 


r = a , 0 = — , 

r 

for some constants a,/3 > 0. The following choice of r and 0 achieves this. 


( 2 . 22 ) 
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Proposition 2.2. Suppose V ,V~^ < 0 and V > 1. Let 


r{x,y) = 1 

lay- 

{a — 

to 

1 

XV- ^ 
v+ ’ 

X 

< 

0 

(2.23) 

[ay- 

{a — 

l)V2/' + 

2x , 

X 

> 

0 

1 


1 + 

a 

y 

). 

X 

< 

0 


9{x,y) = \ 




/ 



1 

(2.24) 

1 

1 + 

OL \ 

—S'-1 1 

J 

+ 27r , 

X 

> 

0 



where tl and tr are given by h2.19\) and and 


a = 


1 + 


1 

V-^ 


/S 


27ra 
ln(/i) ’ 


(2.25) 


where pi is given by \2.1(]\) . Then ^2.23i) and 1(2.241 ) define a continuous bijection from {x,y) G 
to (r, 9) G (0, 0) U (0, oo) x [0, 27r) that satisfies !i2.21\) . Moreover, under this coordinate change 
the restriction of 112.1\} to A is given by 112.22\} with 112.25\} . 


A constructive proof of Proposition 12.21 using (12.6p and (12.71) is given in Appendix lAl Fig. 0] 
illustrates contours of r{x,y) and 9{x,y) as defined in fl2.23l) and fl2.24l) . 

Given fl2.24l) . the negative y-axis corresponds to 6^ = ■ Also by fl2.25p we have 

0 < a < |, with a —)■ 0 as V~V~^ —)■ 1 and a —| as V~V^ —)■ oo (assuming ^ is hxed as the 
limit is taken). Graphically we have found also that 0 < /? < tt. 

The right-hand-side of (I2.22p is dehned everywhere except at r = 0 (the two-fold) and can 
be solved explicitly. Taking r —)■ 0 as t —)■ to as an initial condition, for all t > to the solution to 
fl2.22p is given by 


r{t) = a{t — to) , 9{t) = ln(t — to) + mod 27r , (2.26) 

for some C G [0,27r). Therefore the Filippov solution V’a(t), dehned in 1 12.11 may be written in 
polar coordinates as 

f^a{t) = ir{t),9{t)) = (at, - 

\ a 

2.3 Viable Filippov solutions 

The following result shows that orbits emanating from the two-fold either remain in the repelling 
sliding region R for some time, or they leave R and evolve on A for all later times. 

Proposition 2.3. Let ip{t) = ip{t; 0, 0, 0, 0) be a Filippov solution of ^2.1]) located at the two-fold 
at t = 0. Then either: (i) ip{t) G A for all t > 0, or (ii) ip{t) G R for some t > 0 and ip{t) ^ A 
for allt>0. 


In ( - 


mod 27r 


(2.27) 
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Proof. First, suppose ip{t) G R for some t > 0. Since Filippov solutions to fl2.ip can only enter 
R from the two-fold, either ip{t) G R for all f > 0, in which case clearly (p{t) ^ A for all f > 0, 
or ip{t) E R on time interval (0,6], and at t = h, ip{t) is ejected into either a; < 0 or x > 0. 
Suppose without loss of generality that (p{t) enters x < 0. Then subsequent crossing motion 
(for all f > 0, it can be shown) is described by the return map fl2.8l) starting from (p{b). Since 
(p{b) ^ (, iterations under fl2.8l) do not lie on hence (p(t) ^ A for all f > 0. 

Second, suppose (p(t) ^ R for all f > 0. By solving backward in time, the only way that 
solutions can reach the two-fold without intersecting R is through intersections with (. This is 
because orbits cannot reach the two-fold by evolving backward in time on the attracting sliding 
region. Nor can solutions reach the two-fold by evolving purely in either the left half-space or the 
right half-space because such evolution follows fl2.6l) and fl2.7l) . The only remaining possibility is 
to reach the two-fold via crossing motion, and in view of the saddle-type nature of {y, z) = (0, 0) 
for fl2.8p . this can only be achieved along the unstable eigenvector. Therefore (p{t) G A for all 
t>0. □ 

Now we introduce the notion of a “viable” Filippov solution to represent an orbit that is 
robust for the purpose of forward time numerical simulation. Any orbit that evolves from a point 
in R, simulated using discretisation or by modelling the switch using hysteresis, time-delay, or 
noise, is likely to be immediately ejected into either x < 0 or x > 0. In view of Proposition 
12.31 any simulated orbit that arrives at the two-fold is likely to subsequently evolve on A (or at 
least very near A), as observed in [12]. Indeed, perturbations of a planar two-fold by hysteresis, 
time-delay and noise studied in [18] conhrm that perturbed systems follow Filippov solutions 
provided they do not lie on a repelling sliding region. 

Definition 2.1. A Filippov solution is said to viable if it does not intersect (or travel along) a 




X X 

Figure 4: Contours r{x,y) = constant fl2.23p (panel A), and 0{x,y) = constant fl2.24p (panel 
B). Contours of 6{x,y) are parabolas of the form x = Ky"^, for different values of K. Contours 
of r{x,y) are parabolas of the form x = Ci{y — + C 2 K^, for different values of K, where Ci 

and C 2 are constants that depend on the values of V~ and V~^. In each plot one orbit of fl2.ip is 
shown spiralling away from the origin (shaded blue - colour online). These illustrations use the 
values V~ = —0.5 and V~^ = —2.5. 
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repelling sliding region in forward time. 


For the systems studied in this paper, ip{t;XQ,to) is viable if and only if, for all t > to, 

(p{t-, Xo,to) ^ R- The idea is that viable Filippov solutions are the most relevant in real systems 

or in simulations. The last result of this section follows immediately from Proposition 12.31 and 
provides us with a complete characterisation of viable Filippov solutions to fl2.ip that pass through 
the two-fold. 

Corollary 2.4. Let ip{t) = (p{t; 0, 0, 0, 0) be a Filippov solution of H2.1\) that is located at the 

two-fold at t = 0. Then ip(t) is viable if and only if ^{t) = 'fait), for some a. 

3 Deterministic global dynamics 

The vector held near a generic invisible two-fold in a three-dimensional system can be transformed 
to 

X = /° (3 1) 

\Fs(Jf) + GR(X), x>0 

where Fi and Fr are the constituents of the normal form fl2.2p . and 

aLM(X) ^ {o (x,\y,z\^) ,0{\X\),0(\X\)) , (3.2) 

represent higher order terms. 

As described in [16], there exists a surface A (that coincides with A in the limit |X| —)■ 0) 
on which viable Filippov solutions to fl3.1l) evolve. This surface intersects x = 0 with ?/ > 0 on 
a curve ( that matches ( to hrst order. We let 'fait) denote a viable Filippov solution to fl3.ll) 
with fa{0) = 0 and fa{a) G (. 

Next, in §3.11 we dehne the asymptotic phase for orbits of (13.11) . study the associated isochrons 
in §3.21 and hnd times of return to the discontinuity surface in §3.31 

3.1 The asymptotic phase of an orbit 

Let us suppose that f|3.ip has a stable periodic orbit F of period r such that orbits emanating 
from the two-fold on A approach F. Examples are given below and also in [TS]. In the examples 
considered below, F has exactly two intersections with the discontinuity surface, one with y > 0 
and one with y < 0. We use the intersection point with ?/ > 0 as a reference point with zero 
phase. 

Let (p{t;Xo,to) be a Filippov solution to (13.ip that limits to F as f —)■ cxd. The phase of f, 
relative to a reference time T assumed to be sufficiently large that the orbit lies close to F, will 
be dehned as 


fr = —^- — mod 271 , st = n^^ [x{t) = 0, y{t) > 0] . (3.3) 

Note that sr < T is the previous time at which the orbit lies on x = 0 with |/ > 0. The “mod 27r” 
ensures = [0,27r), but is almost redundant because the orbit is close to F and so T — st is 
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Figure 5: Five isochrons on the surface A for the system (13.Ih with (13. 5 h and (13.dh . 


unlikely to be greater than r. Assuming the forward orbit of Xq is unique and converges to F, 
the asymptotic phase of Xq is dehned as 


(p = lim (pT ■ 

T^oo 


(3.4) 


3.2 Isochrons 

An isochron is a set of points with the same asymptotic phase [T^[20] . Fig. [5]shows hve isochrons 
on A, produced for fl3.1l) with 


V- = -0.5 , V+ = -2.5 , (3.5) 

and 

Gl(X) = Gr{X) = -X . (3.6) 

A was computed by htting a mesh to 1000 numerically computed forward orbits. The stable 
periodic orbit F forms the boundary of this surface. The isochrons were computed by interpolating 
between computed points on the forward orbits and correspond to 0 = ^ for A; = 0,..., 4. More 
sophisticated methods for computing isochrons are discussed in [211 [22]. Each isochron emanates 
transversally from F (as is the case for a generic stable periodic orbit of a smooth system). 

3.3 Times of return to the discontinuity surface 

Here we dehne a return time function / for the curve C,. First note that 'ipa{o-) ^ C by dehnition. 
Then for any a > 0, let t = /(a) > a be the next time at which 'ipait) £ C- Fig. [6] shows a plot of 
/ using (13. 5 p and (13.bh and was computed by numerical simulation. 
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Figure 6: The return time function / for the system (13.11) with (I3.5p and (I3.6p . 


For small a, /(a) ~ /ra, because near the two-fold the higher order terms Gl and Gr have 
little effect and the return time is /la for the normal form fl2.ip . as stated in 1 12.1I For large a, 
ilJaif) is located near F and so /(a) ~ a -|- r. 

In the next section we use / to extrapolate the probability density function for 0^ from points 
near the two-fold to points near F. 


4 Stochastic dynamics and phase randomisation 


In this section we study a stochastic perturbation of (13. ip given by the three-dimensional stochas¬ 
tic differential equation 


dX{t) = 


FL{X{t)) + GL{X{t)) , x(t)<0 

FR{X{t)) + GR{X{t)) , x{t)>Q 


I dt -|- sD dW (t) , 


(4.1) 


where W{t) is a standard three-dimensional Brownian motion, H is a 3 x 3 matrix of constants, 
and 0 < £ 1 represents the noise amplitude. Given a sample solution to fl4.ip . ipe{t;Xo,0), 

and a time T, we define st and 0r by (13. 4 p in the same manner as for the deterministic system 

cu. 

Here we show that the asymptotic phase is highly randomised for sample solutions to fl4.ip 
that pass close to the two-fold before approaching a stable periodic orbit. We provide numerical 
evidence for this in 1 14.11 then derive theoretical approximations to the phase distribution in 1 14.21 
To illustrate the results we use (13. 5 p and two choices for Gr and Gr, namely (13.6p and 

GLiX) = GRiX) = (-a;^ 0) , (4.2) 


because they provide substantially different phase distributions. With different values of V and 
and different choices for Gr and Gr we have observed similar results. 
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Figure 7: Panel A shows a typical sample solution, both in phase space and as a time series, for 
(14.111 with (13-Sp and (I3.6p using D = I and £ = 0.001. Panel B shows an typical sample solution 
using instead (14.211 . Both solutions start at X = (0,1,1), pass near the two-fold X = (0,0,0), 
then approach a stable periodic orbit P. The solutions were computed using the Euler-Maruyama 
method with a step-size of At = 10“^. 


4.1 Sample solutions and phase randomisation 

A sample solution to (14.ip using (I3.5p - (l3.6p is shown in Fig. O-A. The initial point is X = (0,1,1), 
and we denote this solution as ipe(t-, 0,1,1, 0). The deterministic solution, (p(f; 0,1,1, 0), slides into 
the two-fold in a time to = 3.0445 (to four decimal places). The perturbed solution (feit; 0,1,1, 0) 
initially follows a random path close to (p(t; 0,1,1, 0) along the discontinuity surface, and att = to 
is located near the two-fold. While t ^ to, the noise affects <^e(f) qualitatively. After this time 
spirals outward. The initial portion of this spiralling motion is close to the two-fold, and so 
the higher order terms Gl and Gr have little influence. At later times (pe{t) is located relatively 
far from the two-fold. Here Gl and Gr are important and (felt) approaches the stable periodic 
orbit of (13.11) . F. 

We computed 10^ sample solutions analogous to the solution shown in Fig. [7]-A (i.e. with the 
same initial point and parameter values). Fig. [T]-A shows a histogram of the phase 0^ of these 
solutions. Here we used T = 15 because transient dynamics appears to have decayed by this 
time. We notice that the phase is roughly uniformly distributed on [0,27r). 

Fig. [7]-B shows a sample solution using (14.2p instead of (13.61) . As in panel A, the solution 
approaches the two-fold (the deterministic sliding time is to = 2.9763, to four decimal places). 
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then spirals out toward a stable periodic orbit F. Fig. [T]-B shows a histogram of the phase of 
10^ sample solutions, using T = 40. Again the phase is highly random, but in this case the phase 
distribution is not well-approximated by a uniform distribution. A phase oi (pT ^ ^ appears to 
be about twice as likely as a phase of ~ In further simulations (not shown) with different 
values of £ and T we observed similar phase distributions to those shown in Fig. [1] 

In the next section we combine the polar coordinatisation of the local dynamics given in 1 12.21 
with the global return time function / for orbits spiralling out from the two-fold given in 1 13.31 
in order to construct approximations for the phase distributions of Fig. [H 


4.2 The probability distribution of the asymptotic phase 

Let Xq G be such that the deterministic forward orbit of this point, (p{t;Xo,0), is located at 
the two-fold at some time to > 0. For sample solutions (ps{t] Xq, 0), we are only interested in the 
phase (j)T at large values of T (in order to adequately approximate the asymptotic phase 0) but 
our analysis requires considering all T > Iq. 

For any T > Iq, let Pt{st) denote the probability density function for the value of st (the 
previous time of intersection with the discontinuity surface at a point with y > 0). We begin by 
explaining why, for intermediate values of T — to, specihcally e T — fo "C 1, it is suitable to 
assume st — to has a reciprocal probability distribution, that is 


PtM 


1 

■" ( f-HTi-h ) - k) 


(4.3) 


Since the noise amplitude £ is small, if T —to "C 1 then an arbitrary sample solution (pdi) will 
he near the two-fold with high probability. Thus for the purposes of computing st we can ignore 
the higher order terms Gl and Gr. Also if £ -C T — to, then with high probability an arbitrary 
sample solution will he sufficiently far from the two-fold that for the purposes of computing st 
we can ignore the noise. Therefore, here it is suitable to restrict our attention to the normal form 
(EU, and we work with this system in polar coordinates fl2.22p . 

Any solution to 02.221) that limits to the two-fold as f —>■ to (with t > to) is given by 02.26P 
for some G G [0,27r). The system 02.221) is independent of 0 which implies we should treat G as 
a random variable with a uniform distribution on [0, 27r). 

Given G, we can calculate st by using 02.26P to solve 0(t) = 0 for st- We write ^ ln(T —to) = 
27rn -|- 0, for some n G Z and 0 G [0, 27r). By 02.261) . if 0 < G < 27r — 0 then 

— ln(sT — to) -l- G = 27m , (4.4) 

a 

and if 27r — 0 < G < 27r then 

— ln(sT — to) -|- G = 27r(n -|- 1) . (4.5) 

a 

Since G is uniformly distributed, by 04.4p and 04.5p . ln(sT — to) is also uniformly distributed. 
Therefore, in this approximation, st — to has a reciprocal probability distribution as given by 

(lOD . 
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Figure 8: Plots of /”, using n = 10, where / is the return time function dehned in §3.31 Panel A 
corresponds to (13.Oh . for which to ~ 3.0445, and we used T = 15. Panel B corresponds to (14.2p . 
for which tg ~ 2.9763, and we used T = 40. 


As sample solutions (pe{t) continue to evolve outward from the two-fold, since there are no 
other singular points for solutions to encounter, the noise only has the effect of diffusing inter¬ 
section times St by a small amount. Therefore in order to approximate pt for larger values of 
T > to, we can again ignore the noise (as long as T is not so large that the contribution of the 
accumulated small diffusive effects is signihcant). In this approximation, px can be expressed in 
terms of Pf-^(T) by iterating the density under /, specihcally, 

Pt{st) = Pf-^T) (/“^(sr)) ■ (4.6) 

Iterating n times gives 

Pt{st) = Pf-^{T) (/~"(sr)) ■ (4.7) 

We can therefore approximate pt for a large value of T by using fl4.7p with a value of n that is 
sufficiently large that f~'^{T) —to is small, and so the reciprocal probability density function fl4.3p 
can be used for pj-n^T)- This is the manner by which we obtained the two approximations shown 
in Fig. m using n = 10 in both cases. The iterative procedure fl4.7l) was performed numerically 
because analytic expressions for / appear to be unavailable. 

Fig. [8] shows plots of /"■ for both fl3.6p and fl4.2p . We note that if /""(a) is an affine function 
of ln(a) then our procedure generates a uniform distribution for pt- Indeed in Fig. |8]-A /"■ is 
indistinguishable from an affine function on the given scale and so in Fig. [I]- A the distribution 
is approximately uniform. In Fig. [SJ-B, /” has a noticeable nonlinearity and for this reason the 
distribution in Fig. [T]-B is signihcantly non-uniform. 

5 Desynchronising a collection of smooth oscillators 

Suppose we wish to break the synchrony of a collection of coordinated oscillators. In [7] this 
is achieved by applying a control that pushes the state of the oscillators toward an unstable 
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equilibrium, where small random perturbations efficiently randomise the phase of the oscillators. 
Upon removal of the control action, each oscillator returns to the original regular periodic motion 
but now the oscillators are desynchronised. Here we propose a “fast” method to achieve this, 
using a two-fold as the phase singularity rather than an equilibrium, in which case there is no 
slowing of the dynamics as the oscillators return to periodic motion. 

To illustrate our method we use the Hopf bifurcation normal form 

(h, y) = F{x, y) = [x-y- x{x^ + y‘^),x + y- y{x^ + y'^)) , (5.1) 

and apply a discontinuous control that creates a two-fold singularity. In principle this can be 
achieved with any system exhibiting a stable periodic orbit, including excitable systems such as 
the FitzHugh-Nagumo equations [7]. 

Consider the two-dimensional stochastic differential equation 

(^dx{t), dy{t)^ = (^F{x{t),y{t)) + H{t — — t)c{t)) dt + e dW(t) , (5.2) 

where W{t) is a standard two-dimensional Brownian motion, and 0 < e < 1 is the noise ampli¬ 
tude. H is the Heaviside function and ti and t 2 are the start and end times of a control action 
given by 

f(ait,a 2 ), a ;<0 
cit) = \(, . ^ n ’ 

[^( 03 ^, 04 ), a; > 0 

where oi, 02 , 03 , 04 G M are control parameters. 

By treating the time t as a third dynamic variable, i.e. with t = 1, the system is three- 
dimensional and is piecewise-smooth when ti < t < t 2 with the discontinuity surface x = 0. The 
fold lines on a; = 0 (the lines where x = 0) are y = ait and y = so the origin (x, y, t) = (0, 0, 0) 
is a two-fold. Via elementary calculations we hnd we need 02 < oi, 03 < 04 and ai 7 ^ 03 in order 
for the two-fold to be generic and for both folds to be invisible, as in fl2.1l) . We also require 
Oi < 03 in order to have V~, V~^ < 0 and V~V~^ > 1 fl2.5p . In summary, we require 


02 < Oi < 03 < 04 . (5.4) 

Fig. [9] illustrates the phase randomisation using 

Oi = —0.2 , 02 = —1 , 0,3 = 0.2 , 04 = 1 , (5.5) 

and 

ti = -5 , t2 = 2.5 . (5.6) 

In simulations we found that with relatively small values of the o,, a large negative value of ti is 
required in order for the control to direct orbits into the two-fold, and that with large values of 
Oj the distributions of the asymptotic phase were substantially non-uniform. 

6 Discussion 

In this paper we have studied orbits that pass through an invisible two-fold (Teixeira singularity) 
then limit to a stable periodic orbit. We considered the presence of small noise in order to 
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Figure 9: Panel A shows five sample solutions to fl5.2p using flS.Sp . 05.61) and £ = 0.001, and the 
initial condition {x,y) = (1,0) at t = —15. The control 05.3p causes the solutions to pass close 
to a two-fold at t = 0, effectively randomising the phase. Panel B shows one of the solutions 
in phase space, with the location of the solution highlighted at key times. The sample solutions 
were computed using the Euler-Maruyama method with a step-size of At = 10“®. 


remove the ambiguity of evolution through the two-fold. We found that the phase of the limiting 
periodic motion is highly randomised and that the probability distribution of the phase depends 
crucially on the nature of the nonlinear dynamics experienced by orbits between the two-fold and 
the periodic orbit. By using polar coordinates to describe the motion of orbits as they initially 
depart from the two-fold, and a one-dimensional return map to describe the global dynamics 
approaching the periodic orbit, we constructed approximations to the probability density function 
of the asymptotic phase that matched well to the results of numerical simulations. Fig. [H 

In 1|5] we showed how a simple discontinuous control law can generate a two-fold in a smooth 
system. We propose that in this fashion the two-fold can be used to desynchronise a collection 
of oscillators and that this may have the following advantages to desynchronisation using an 
unstable eqnilibrium. 

i) With an unstable equilibrium, the control action must direct the state of each oscillator 
to near the equilibrium, and the effectiveness of the randomisation correlates with the 
accuracy of this action. With a two-fold, however, the control does not need not be as 
precise because it only needs to direct orbits to the basin of attraction of the two-fold. 

ii) With an unstable equilibrium, the effect of the randomisation is due to inherent stochas- 
ticity in the system and any artificial randomness in the control law. With a two-fold, the 
randomisation is inherent in the ambiguity of forward evolution through the two-fold. 

hi) Orbits dwell near unstable equilibria and so can take a relatively long time to retnrn to 
regular periodic motion. With a two-fold there is no slowing of the dynamics, hence we 
refer to our desynchronisation as “fast” phase randomisation. Indeed, in Fig. [9] we see 
that solutions return quickly to approximately regular periodic motion after the control is 
turned off. Whether this truly constitutes fast desynchronisation requires a study of the 
full time the control action is required to act for. This, and a study of practical control 
actions, remain for future work. 
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Teixeira’s paper [9] inspired a legacy of intrigue around a specific case of the two-fold singu¬ 
larity, first studied more generally in [ 8 ]. Interest has only increased as its role as a determinacy- 
breaking singularity has become more clear [TUlfTB] . This paper begins to reveal the practical side 
of these insights, including how two-folds may manifest in physical systems, and how they might 
be put to use as tools for control. Little is still known about where two-folds appear naturally 
in physical systems. Similarly little is known about what applications they may have in control 
systems, but a role as a phase randomiser is suggested here. Perhaps the obvious next step is 
to design implementable control circuits to investigate the practical obstacles and opportunities 
that they present. 


A Proof of Proposition 12.2 


Proof. The ?/-value of '0a(a) is given by (I2.16p . and so by (I2.27p the r-value of 'ipaia) is equal to 
aa. Since 'ipaict) G C) the y and r-values of V’a(a) are the same, see fl2.2ip . hence a is given by 
(1125]). 

At the start of 1 12.21 we showed that, given x < 0 and ?/ G M, if > 0 is given by (12.191) and 
Vo = y — V~tl, then evolving (12.11) forward from (0,|/o,7|/o) for a time takes us to the point 
{x,y,z), where x = ^'P{y,z) (I2.14p . Moreover, this dehnes a bijection from the region x < 0, 
y eM to the region j/o > 0, 0 < tl < — 277 / 0 - 

By (I2.22p . evolving (r, 9) = (7/0, 0) for a time tl takes us to the point 

(^r,e) = ^yo + axL, + 1^^ . (A.l) 

This is a bijection from 7/0 > 0, 0 < < — 277/0 tor>0, 0<6'<^ ln(l — 2 q! 7) (assuming 

/3 > 0). By substituting (I2.19P and 7/0 = 7 / — V~tl into (lA.ip and simplifying, we arrive at (I2.23P 
and (I2.24P for x < 0. 

Similarly, given x > 0 and 7 / G M, let < 0 be given by (I2.20|) and yo = y — tr. This is a 
bijection to 7/0 > 0, —2yo < tr < 0, and evolving (12.11) backward from ( 0 , 7 / 0 , 77 / 0 ) for a time tr 
takes us to (x, 7 /, z), where x = r^E{y, z) (I2.14p . 

In polar coordinates evolving (r, 9) = (7/0, 0) for a time tr takes us to 


{r,9) 


7/0 axR, - In 
a 




(A. 2 ) 


which is a bijection to r > 0, - ln(l — 2a) + 27: < 9 < 2 tt. Substituting (I2.20p and yo = y — tr 
into (IA.2p leads to (I2.23P and (I2.24p for x > 0. 

Finally, by matching the limiting left and right values for 9 on the negative //-axis as given 
by (I2.24p . we obtain 

9 = — ln(l — 207 ) = — ln(l — 2a) -|- 27r . (A.3) 

a a 


m 


By using the identity y = by solving (IA.3P for (3 we obtain the expression for (3 given i 

mB- 

This completes the proof because the last calculation provides continuity of 9{x,y), and by 
construction r and 9 satisfy the ordinary differential equations (I2.22p . □ 
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